##########################################
########################################## 
## STUDY 6
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

# load data for study
data_raw <- haven::read_dta(here("Data", "study6.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df <- data_raw %>% 
  mutate(female = case_when(
    q33 == 2 ~ 1,
    q33 == 1 ~ 0,
    TRUE ~ NA_real_
  ),
  
  #Behavior in Hawk-Dove game
  hawk01 = case_when(
    q44 == 1 ~ 0,
    q44 == 2 ~ 1,
    TRUE ~ NA_real_
  )) %>% 
  rename(age = q34,
         race = q35,
         edu = q36) %>% 
  
  #Party Identification (1 = Democrat; 0 = Republican)
  mutate(repdem = case_when(
    q38_1 == 1 | q38_1 == 2 | q38_3 == 1 ~ 0,
    q38_2 == 1 | q38_2 == 2 | q38_3 == 2 ~ 1,
    TRUE ~ NA_real_)) %>% 
  
  #Dummy indicating whether participant is Independent or not
  mutate(independents = case_when(
    q38_3 == 3 ~ 1,
    TRUE ~ NA_real_)) %>% 
  
  #Need for Chaos
  mutate_at(vars(need_chaos_1:need_chaos_7), as.numeric) %>% 
  mutate(need_chaos = rowMeans(select(., c("need_chaos_1", "need_chaos_2", "need_chaos_3", 
                                           "need_chaos_4", "need_chaos_5", "need_chaos_6", "need_chaos_7"))), 
         na.rm = TRUE) %>%
  #Loneliness
  mutate_at(vars(lonely_1:lonely_7), as.numeric) %>% 
  mutate(lonely_5 = flip(lonely_5)) %>% 
  mutate(loneliness = rowMeans(select(., c("lonely_1", "lonely_2", "lonely_3", 
                                           "lonely_4", "lonely_5", "lonely_6", "lonely_7"))), 
         na.rm = TRUE) %>% 
  
  #Social Dominance Orientation
  mutate_at(vars(sdo_1:sdo_8), as.numeric) %>% 
  mutate_at(vars(sdo_3, sdo_4, sdo_7, sdo_8), ~flip(.)) %>% 
  mutate(SDO = rowMeans(select(., sdo_1:sdo_8), 
                        na.rm = TRUE)) %>% 
  
  #Psychopathy
  mutate_at(vars(triad1_1:triad3_9), as.numeric) %>% 
  mutate_at(vars(triad3_2, triad3_7), ~flip(.)) %>%  
  mutate(dark_triad_3 = rowMeans(select(.,triad3_1:triad3_9)), 
         na.rm = TRUE) %>% 
  
  #Status-Driven Risk-Taking
  mutate_at(vars(sdrt_1:sdrt_14), as.numeric) %>% 
  mutate_at(vars(sdrt_1, sdrt_7, sdrt_8, sdrt_12, sdrt_14), ~flip(.)) %>% 
  mutate(sdrt = rowMeans(select(.,sdrt_1:sdrt_14), 
                         na.rm = TRUE)) %>% 
  
  #Competitive Jungle Worldview
  mutate_at(vars(comp_jungle_1:comp_jungle_12), as.numeric) %>%   
  mutate_at(vars(comp_jungle_2, comp_jungle_4, comp_jungle_6, comp_jungle_8, comp_jungle_11, comp_jungle_12), ~flip(.)) %>%  
  mutate(comp_jungle = rowMeans(select(.,comp_jungle_1:comp_jungle_12)), 
         na.rm = TRUE)  %>% 
  
  #Social Ladder Placement
  mutate_at(vars(ladder), as.numeric) %>% 
  
  #Political Rumors: Believing and Sharing of DEMocratic and REPublican rumors
  mutate(true_dem = rowMeans(select(.,c("dem1_1", "dem2_1", "dem3_1"))), 
         na.rm = TRUE) %>% 
  mutate(share_dem = rowMeans(select(.,c("dem1_2", "dem2_2", "dem3_2"))), 
         na.rm = TRUE) %>% 
  mutate(true_rep = rowMeans(select(.,c("rep1_1", "rep2_1", "rep3_1"))), 
         na.rm = TRUE) %>% 
  mutate(share_rep = rowMeans(select(.,c("rep1_2", "rep2_2", "rep3_2"))), 
         na.rm = TRUE) %>% 
  #Standardize variables 
  mutate(across(.cols = c(female, age, race, edu, ladder ,need_chaos, loneliness, SDO, dark_triad_3, sdrt, comp_jungle, true_dem, share_dem, true_rep, share_rep),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_z_score")) %>% 
  
  #Scale variables to range from 0 to 1
  mutate(across(.cols = c(female, age, race, edu, ladder, need_chaos, loneliness, SDO, dark_triad_3, sdrt, comp_jungle, true_dem, share_dem, true_rep, share_rep),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}01"))

##########################################
########################################## 
## Sample Description 
## (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
#Gender
print(paste(round(mean(df$female , na.rm =TRUE), 2)*100 , "percent of participants are female"))
#Age
print(paste("Mean age of participants:" , round(mean(df$age , na.rm =TRUE), 2), 
            ", SD:", round(sd(df$age, na.rm = TRUE) , 2) ))
#Education
print(paste("Less than high school:", round(prop.table(table(df$edu)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df$edu)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df$edu)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df$edu)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df$edu)), 2)[5]*100, "percent"))
print(paste("Graduate or professional degree:", round(prop.table(table(df$edu)), 2)[6]*100, "percent"))
#Ethnicity
print(paste("White/caucasian:", round(prop.table(table(df$race)), 2)[2]*100, "percent"))

##########################################
########################################## 
## Factor Analysis, etc.
## ## (See Supplementary Materials SM2)
##########################################
########################################## 
#Cronbach's alpha: Need for Chaos
nfc <- df %>% 
  select(need_chaos_1:need_chaos_7) %>% 
  psych::alpha(. , check.keys = TRUE) 
print(paste("Cronbach's alpha, Need for Chaos:" , round(nfc$total$raw_alpha,2)))  

#Confirmatory factor analysis: Need for Chaos
cfa <- df %>% 
  select(need_chaos_1:need_chaos_7) %>% 
  na.omit()

#One-Factor Model (Table S7)
m1  <- ' nfc_oneDim  =~ need_chaos_1+need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

#Cronbach's alpha, mean, standard deviation: Loneliness
alone <- df %>% 
  select(lonely_1:lonely_7) %>% 
  psych::alpha(. , check.keys = TRUE) 
print(paste("Cronbach's alpha, Loneliness:" , round(alone$total$raw_alpha,2)))  

mean(df$loneliness01, na.rm = T)
sd(df$loneliness01, na.rm = T)

#Mean & standard deviation: Social Ladder Placement
mean(df$ladder01, na.rm = T)
sd(df$ladder01, na.rm = T)

#Cronbach's alpha, mean, standard deviation: Status-Driven Risk-Taking
status_risk <- df %>% 
  select(sdrt_1:sdrt_14) %>% 
  psych::alpha(. , check.keys = TRUE) 
print(paste("Cronbach's alpha, Loneliness:" , round(status_risk$total$raw_alpha,2)))  

mean(df$sdrt01, na.rm = T)
sd(df$sdrt01, na.rm = T)

##########################################
########################################## 
## Remove 5% with highest Need for Chaos? (SM8)
## When set to TRUE, used for producing Fig S8, Fig S9, Fig S12
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df <- df %>% 
    filter(need_chaos_z_score < quantile(need_chaos_z_score, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 3 of Main Text
##########################################
########################################## 
m_share_dem <- lm(share_dem01 ~ need_chaos_z_score + repdem  +  female + edu_z_score + age_z_score + race , data = df )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_z_score" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study6")

m_share_rep <- lm(share_rep01 ~ need_chaos_z_score + repdem  +  female + edu_z_score + age_z_score + race , data = df )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_z_score" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study6")

m_true_dem <- lm(true_dem01 ~ need_chaos_z_score + repdem  +  female + edu_z_score + age_z_score + race , data = df )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_z_score" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study6")

m_true_rep <- lm(true_rep01 ~ need_chaos_z_score + repdem  +  female + edu_z_score + age_z_score + race , data = df )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_z_score" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study6")

study6_coef <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study6_coef,file=here("Data_recoded","study6_fig3.R"))

#regression table for supplemental materials (Table S13)
stargazer(m_share_dem , m_share_rep , m_true_dem , m_true_rep , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Share - Dem", "Share - Rep" , "True - Dem" , "True - Rep"),
          covariate.labels = c("Need for Chaos (std.)" , "Party ID (1 = Dem.)" , "Female" , "Education" , "Age  (std.)" , "Ethnicity (1 = non-white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

##########################################
########################################## 
## Analysis - Fig. 4 of Main Text
##########################################
########################################## 
#Fig 4: Democratic Voters 
df_dem <- df %>% 
  filter(repdem == 1 )

m_share_dem <- lm(share_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_dem )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study6")

m_share_rep <- lm(share_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_dem )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study6")

m_true_dem <- lm(true_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_dem )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study6")

m_true_rep <- lm(true_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_dem )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study6")

study6_coef_dem <- rbind(share_dem, share_rep, true_dem, true_rep)

#Fig 4: Republican Voters 
df_rep <- df %>% 
  filter(repdem == 0 )

m_share_dem <- lm(share_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_rep )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study6")

m_share_rep <- lm(share_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_rep )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study6")

m_true_dem <- lm(true_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_rep )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study6")

m_true_rep <- lm(true_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_rep )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study6")

study6_coef_rep <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study6_coef_dem,file=here("Data_recoded","study6_fig4_dem.R"))
save(study6_coef_rep,file=here("Data_recoded","study6_fig4_rep.R"))

##########################################
########################################## 
## Analysis - Fig. S2 of Supplementary Materials
##########################################
########################################## 
df_independents <- df %>% 
  filter(independents == 1 )

m_share_dem <- lm(share_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_independents )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study6")

m_share_rep <- lm(share_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_independents )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study6")

m_true_dem <- lm(true_dem01 ~ need_chaos_z_score +   female + edu_z_score + age_z_score + race , data = df_independents )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study6")

m_true_rep <- lm(true_rep01 ~ need_chaos_z_score +  female + edu_z_score + age_z_score + race , data = df_independents )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study6")

study6_coef_independents <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study6_coef_independents,file=here("Data_recoded","study6_figS2.R"))

##########################################
########################################## 
## Analysis - Supplementary Materials SM2
## Association between Need for Chaos and
## other measures of anti-social traits
##########################################
########################################## 

###############################
## Need for Chaos versus Psychopathy
###############################
cfa <- df %>% 
  select(need_chaos_1:need_chaos_7 , triad3_1: triad3_9 ) %>% 
  na.omit()

#Confirmatory factor analysis: One dimension
m1  <- ' nfc_oneDim  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7+triad3_1+triad3_2+triad3_3+triad3_4+triad3_5+triad3_6+triad3_7+triad3_8+triad3_9
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5
        triad3_2~~triad3_7
        triad3_5~~triad3_6'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

#Confirmatory factor analysis: Two dimensions
m2  <- ' chaos  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5        
        psych  =~ triad3_1 +triad3_2+triad3_3+triad3_4+triad3_5+triad3_6+triad3_7+triad3_8+triad3_9
        triad3_2~~triad3_7
        triad3_5~~triad3_6'
m2 <- lavaan::cfa(m2, data=cfa) 
m2_fit <- summary(m2 , fit.measures=TRUE) 

print(paste("Two-factor model, Chi-square:" , round(m2_fit$FIT["chisq"],2 )))
print(paste("Two-factor model, CFI:" , round(m2_fit$FIT["cfi"],2 )))
print(paste("Two-factor model, TLI:" , round(m2_fit$FIT["tli"],2 )))
print(paste("Two-factor model, RMSEA:" , round(m2_fit$FIT["rmsea"],2 )))

#Compare models
compare <- anova(m1, m2, method = "satorra.bentler.2010")

print(paste("Chi-square difference:" , round(compare[2,5],2)))
print(paste("Chi-square difference, p-value:" , compare[2,7]))

###############################
## Need for Chaos versus Social Dominance Orientation
###############################
cfa <- df %>% 
  select(need_chaos_1:need_chaos_7, sdo_1: sdo_8) %>% 
  na.omit()

#Confirmatory factor analysis: One dimension
m1  <- ' nfc_oneDim  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7 + sdo_1+ sdo_2+ sdo_3+ sdo_4+ sdo_5+ sdo_6+ sdo_7+ sdo_8
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5
        sdo_4~~sdo_7
        sdo_4~~sdo_8
        sdo_7~~sdo_8'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

#Confirmatory factor analysis: Two dimensions
m2  <- ' chaos  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5        
        sdo  =~ sdo_1+ sdo_2+ sdo_3+ sdo_4+ sdo_5+ sdo_6+ sdo_7+ sdo_8
        sdo_4~~sdo_7
        sdo_4~~sdo_8
        sdo_7~~sdo_8'
m2 <- lavaan::cfa(m2, data=cfa) 
m2_fit <- summary(m2 , fit.measures=TRUE) 

print(paste("Two-factor model, Chi-square:" , round(m2_fit$FIT["chisq"],2 )))
print(paste("Two-factor model, CFI:" , round(m2_fit$FIT["cfi"],2 )))
print(paste("Two-factor model, TLI:" , round(m2_fit$FIT["tli"],2 )))
print(paste("Two-factor model, RMSEA:" , round(m2_fit$FIT["rmsea"],2 )))

#Compare models
compare <- anova(m1, m2, method = "satorra.bentler.2010")

print(paste("Chi-square difference:" , round(compare[2,5],2)))
print(paste("Chi-square difference, p-value:" , compare[2,7]))

###############################
## Need for Chaos versus Competitive Jungle Worldview
###############################
cfa <- df %>% 
  select(need_chaos_1:need_chaos_7, comp_jungle_1: comp_jungle_12) %>% 
  na.omit()

#Confirmatory factor analysis: One dimension
m1  <- ' nfc_oneDim  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7 +comp_jungle_1+comp_jungle_2+comp_jungle_3+comp_jungle_4+comp_jungle_5+comp_jungle_6+comp_jungle_7+comp_jungle_8+comp_jungle_9+comp_jungle_10+comp_jungle_11+comp_jungle_12
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

#Confirmatory factor analysis: Two dimensions
m2  <- ' chaos  =~ need_chaos_1 +need_chaos_2+need_chaos_3+need_chaos_4+need_chaos_5+need_chaos_6+need_chaos_7
        need_chaos_3~~need_chaos_4
        need_chaos_3~~need_chaos_5
        need_chaos_4~~need_chaos_5        
        cjw  =~ comp_jungle_1+comp_jungle_2+comp_jungle_3+comp_jungle_4+comp_jungle_5+comp_jungle_6+comp_jungle_7+comp_jungle_8+comp_jungle_9+comp_jungle_10+comp_jungle_11+comp_jungle_12
       '
m2 <- lavaan::cfa(m2, data=cfa) 
m2_fit <- summary(m2 , fit.measures=TRUE) 

print(paste("Two-factor model, Chi-square:" , round(m2_fit$FIT["chisq"],2 )))
print(paste("Two-factor model, CFI:" , round(m2_fit$FIT["cfi"],2 )))
print(paste("Two-factor model, TLI:" , round(m2_fit$FIT["tli"],2 )))
print(paste("Two-factor model, RMSEA:" , round(m2_fit$FIT["rmsea"],2 )))

#Compare models
compare <- anova(m1, m2, method = "satorra.bentler.2010")

print(paste("Chi-square difference:" , round(compare[2,5],2)))
print(paste("Chi-square difference, p-value:" , compare[2,7]))

##########################################
########################################## 
## Analysis - Supplementary Materials SM8
## Association between Need for Chaos and
## behavior in Hawk-Dove game
##########################################
########################################## 
df_reduced_hawk_dove <- df %>% 
  select(need_chaos01, need_chaos_z_score, repdem, female, edu_z_score, age_z_score, race, hawk01) %>% 
  na.omit()

reg_hawk_dove <- lm(hawk01 ~ need_chaos_z_score + female + edu_z_score + age_z_score + race,
                    data = df_reduced_hawk_dove)

#Table S30
stargazer::stargazer(reg_hawk_dove,
                     align=TRUE,
                     no.space=TRUE,
                     dep.var.labels.include = FALSE,
                     keep.stat = c("n","ser"),
                     keep=c("need_chaos_z_score", "female", "edu_z_score",  "age_z_score", "race", "Constant"),
                     covariate.labels = c("Need for Chaos (std.)", "Gender" , "Education", "Age", "Ethnicity",
                                          "Intercept"))

##########################################
########################################## 
## Analysis - Fig. 7 of Main Text
##########################################
########################################## 
df_reduced <- df %>% 
  select(need_chaos01, need_chaos_z_score, repdem, female, edu_z_score, age_z_score, race, sdrt_z_score, loneliness_z_score, ladder_z_score) %>% 
  na.omit() %>% data.frame()

##Main effects of social status and marginalization on Need for Chaos
reg1 <- lm(need_chaos01 ~ sdrt_z_score + loneliness_z_score + ladder_z_score + female + edu_z_score + age_z_score,
           data = df_reduced)
reg2 <- lm(need_chaos01 ~ sdrt_z_score + female + edu_z_score + age_z_score,
           data = df_reduced)
reg3 <- lm(need_chaos01 ~ loneliness_z_score + female + edu_z_score + age_z_score,
           data = df_reduced)
reg4 <- lm(need_chaos01 ~  ladder_z_score + female + edu_z_score + age_z_score,
           data = df_reduced)
reg5 <- lm(need_chaos01 ~ sdrt_z_score*loneliness_z_score + sdrt_z_score*ladder_z_score + sdrt_z_score*female + sdrt_z_score*edu_z_score + sdrt_z_score*age_z_score +
             loneliness_z_score*female + loneliness_z_score*edu_z_score + loneliness_z_score*age_z_score +
             ladder_z_score*female + ladder_z_score*edu_z_score + ladder_z_score*age_z_score,
           data = df_reduced)

##Panel A of Fig. 7: Main effects on Need for Chaos
df_plot <- tidy(reg1) %>% 
  filter(term %in% c("sdrt_z_score", "loneliness_z_score", "ladder_z_score"))

df_plot$term <- factor(df_plot$term, levels = df_plot$term[order(df_plot$estimate)])
df_plot$Variable <- plyr::revalue(df_plot$term , c("sdrt_z_score"="Status-Driven Risk-Taking", 
                                                   "loneliness_z_score"="Loneliness", 
                                                   "ladder_z_score" = "Social Ladder"))

text <- paste("beta == ", round(df_plot$estimate , digits = 2) )

fig6_a <- ggplot(df_plot) +
  geom_point(aes(term, estimate, shape = Variable), size = 3) +
  geom_linerange(aes(x = term, y = estimate, ymin = estimate-2*std.error, ymax = estimate+2*std.error) , 
                 size = 2, alpha = .2) +
  ylab("Status and Marginalization on Need for Chaos") +
  ggtitle("") +
  xlab("") +  
  theme(
    legend.background = element_rect(fill = 'white', linetype='solid'),
    legend.justification=c(1,0), legend.position=c(.98,0.04),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.text=element_text(size=11),
    axis.text.x = element_text(size = 11 , color = "black"),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)) + 
  geom_text(aes(x = term, y = estimate),
            label=text, 
            nudge_y = 0.0, nudge_x = .4, 
            check_overlap = T , parse = TRUE, size  = 6
  ) + 
  coord_flip() + geom_hline(yintercept=0) +
  labs(color = "Variable") +
  guides(shape = guide_legend(reverse = T)) + 
  scale_y_continuous(breaks = c(-.1,-.05,0,.05,.1,.15,.2,.3), 
                     labels = c("-.1","-.05","0",".05",".1",".15",".2",".3"),
                     limits = c(-.1,.15))

##Panel B & C of Fig. 7: Interaction effects

#Interaction between Social Ladder Placement and Status-Driven Risk-Taking on Need for Chaos
fig6_b <-interflex(Y = "need_chaos01", X = "ladder_z_score", D="sdrt_z_score", Z = c("female","edu_z_score","age_z_score", "loneliness_z_score"),
                   data = df_reduced, 
                   na.rm = TRUE,
                   nbins = 5,
                   Xdistr = "histogram",
                   estimator = "binning", 
                   vcov.type = "robust", 
                   xlab = "Social Ladder Placement", 
                   ylab = "Status-Driven Risk-Taking on NFC",
                   show.grid = TRUE,
                   theme.bw = TRUE)
fig6_b <- fig6_b$figure

#Interaction between Loneliness and Status-Driven Risk-Taking on Need for Chaos
fig6_c <-interflex(Y = "need_chaos01", X = "loneliness_z_score", D="sdrt_z_score", Z = c("female","edu_z_score","age_z_score", "ladder_z_score"),
                   data = df_reduced, 
                   na.rm = TRUE,
                   nbins = 5,
                   Xdistr = "histogram",
                   estimator = "binning", 
                   vcov.type = "robust", 
                   xlab = "Loneliness", 
                   ylab = "Status-Driven Risk-Taking on NFC",
                   show.grid = TRUE,
                   theme.bw = TRUE)
fig6_c <- fig6_c$figure

# Combine plots
cowplot::ggdraw() +
  draw_plot(fig6_a, x = 0, y = .5, width = 1, height = .46) +
  draw_plot(fig6_c, x = 0.5, y = 0, width = .5, height = .51) +
  draw_plot(fig6_b, x = 0, y = 0, width = .5, height = .51) +
  draw_plot_label(label = c("A", "B", "C"), size = 14,
                  x = c(0, 0, 0.5), y = c(.93, 0.5, 0.5)) 

#Save
ggsave(here::here("figures", "fig7.pdf"),
       width = 9.2, height =8.2)

##Regression table for Supplemental Materials (Table S20)
stargazer::stargazer(reg1,reg5,
                     align=TRUE,
                     no.space=TRUE,
                     dep.var.labels.include = FALSE,
                     keep.stat = c("n","ser"),
                     font.size = "small", 
                     covariate.labels = c("Status-Driven Risk-Taking (SDRT)", "Loneliness" , "Perceived Status",
                                          "Women","Education", "Age",
                                          "SDRT X Loneliness", "SDRT X Perceived Status", "SDRT X Women", "SDRT X Education",
                                          "SDRT X Age", "Loneliness X Women", "Loneliness X Education",
                                          "Loneliness X Age", "Perceived Status X Women", "Perceived Status X Education",
                                          "Perceived Status X Age", "Intercept"))





